cpi uiaii 



Dynamics of gene expression and the regulatory inference problem 



Johannes Berg 



00 
O 
O 



in 



o 



(N 
> 

cn 
o 



Institut fur Theoretische Physik, Universitdt zu Koln 

Zulpicherstr. 77, 50937 Koln, Germany 

and 

Kavli Institute for Theoretical Physics, 
UCSB, Santa Barbara, CA 93106-4030 



PACS 87. 16.Yc - Regulatory genetic and chemical networks 
PACS 87 . 10 . Mn - Stochastic modeling 
PACS 87.16.dj - Dynamics and fluctuations 

Abstract. - From the response to external stimuli to cell division and death, the dynamics of 
living cells is based on the expression of specific genes at specific times. The decision when to 
express a gene is implemented by the binding and unbinding of transcription factor molecules to 
regulatory DNA. Here, we construct stochastic models of gene expression dynamics and test them 
on experimental time-series data of messenger-RNA concentrations. The models are used to infer 
biophysical parameters of gene transcription, including the statistics of transcription factor-DNA 
binding and the target genes controlled by a given transcription factor. 



Introduction. — The primary step in the production 
of a protein is the transcription of a gene from the DNA 
template to a m(essenger)-RNA molecule. Gene transcrip- 
tion is carefully controlled, allowing the cell to respond to 
situations requiring different proteins: for instance, en- 
zymes are produced when nutrients need to be processed 
and repair proteins are assembled to respond to cell dam- 
age. 

The cellular information processing to achieve this con- 
trol has a concrete biophysical basis. Transcription factor 
molecules locate and bind to specific sites on DNA found 
in the regulatory region in front of a gene. The transcrip- 
tion factors then attract or repel the molecular machinery 
that reads off the gene (see Fig. 1(a)). The resulting gene 
product may itself be a transcription factor, enhancing or 
suppressing the transcription of further genes. The dy- 
namics of gene expression is thus governed by a complex 
web of molecular interactions. Changes in the regulatory 
interactions (caused for instance by mutations in a tran- 
scription factor binding site) disrupt the cellular dynamics 
and can lead to serious defects, including cancer. 

Which transcription factor targets which gene can in 
principle be determined experimentally, by investigating 
the physical binding of transcription factors to the regu- 
latory DNA of a target gene. High-throughput methods 
have resulted in lists of potential target genes for many 



transcription factors [1,2], which, however, contain many 
false positives and false negative results. Binding sites for 
transcription factors can also be identified computation- 
ally by searching regulatory DNA for statistically overrep- 
resented subsequences [3,4]. 

Both approaches, however, tell nothing about the func- 
tion of a transcription factor (for instance its activity as an 
inhibitor or an enhancer of gene transcription), or about 
its interplay with other transcription factors. Moreover, 
some transcription factors (called co-factors) do not bind 
directly to DNA but to other transcription factors, and 
thus do not have binding sites on DNA. As a result, even 
in well-studied organisms such as yeast, the function and 
target genes of about half the transcription factors remain 
unknown [5]. 

The development of microarrays has made it possible 
to simultaneously measure the expression levels (mRNA 
concentrations) of all genes of an organism (see Fig. 1(b)). 
This allows to probe the response of the regulatory net- 
work to external perturbations, or to detect correlations 
of expression levels of different genes [7] . 

The effort to deduce the input-output relations of regu- 
latory networks from experimental gene expression data, 
however, is impeded by the large number of such rela- 
tions which are compatible with a given set of expres- 
sion levels [81. For instance, an observed correlation be- 
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Fig. 1: Biophysics of gene transcription and its experi- 
mental measurement, a) Transcription factors are proteins 
which bind to specific sites on DNA. These binding sites are 
located in the so-called regulatory region of a gene, which for 
yeast typically extends over several hundred nucleotides before 
the starting point of transcription (marked with the arrow). 
A fully assembled set of transcription factors will attract and 
activate the machinery (called RNA-polymerase) which tran- 
scribes the DNA template to an mRNA molecule. In a further 
step, the mRNA molecule will be translated to a protein, 
b) DNA-microarrays can simultaneously measure mRNA ex- 
pression levels in a sample for all genes of an organism. A sam- 
ple is taken from a population of cells in a certain state, so even 
for low expression levels of a gene, the number of molecules is 
macroscopic. (Single-cell fluctuations are averaged out in this 
process.) mRNA in a sample is reversely transcribed to DNA 
and deposited on a chip, which has a strand of complementary 
DNA grafted onto its surface. The DNA in the sample will 
bind to its complement on the chip. Fluorescent labelling of 
the DNA in the sample allows to determine the original mRNA 
concentration relative to a reference sample. By photolitho- 
graphic techniques, hundreds of thousands of such DNA-probes 
can be placed on a single microarray-chip. By now over 10 
genome- wide experimental measurements of mRNA-expression 
levels have been performed for different organisms, under dif- 
ferent experimental conditions, or for different tissues [6]. 

tween mRNA concentration of two genes may stem from 
a causal relationship (one gene regulates the expression of 
the other), or from both genes sharing a common regu- 
lator [7]. As a result, many machine learning approaches 
to the network reconstruction problem, such as Baycsian 
probabilistic models (see [9,10] for reviews), or linear re- 
gression [11] are well-known to be afflicted with the curse 
of dimensionality [8, 12]. 

Here we take a biophysics approach to gene transcrip- 
tion. The inference of the biophysical machinery underly- 



ing gene transcription from its output is a challenging in- 
verse problem akin to unravelling the principles of atomic 
physics from spectroscopy or the events in the early uni- 
verse from cosmic radiation. In this paper, the statistics 
of gene expression dynamics is deduced from experimen- 
tal data. Based on this statistics, a Langevin model for 
the non-equilibrium dynamics of mRNA concentrations is 
used to infer basic biophysical parameters of expression 
dynamics, including mRNA degradation constants, the 
statistics of transcription factor-DNA interactions, and 
regulatory interactions. 

Modelling gene expression dynamics. The ba- 
sic processes driving the concentration of mRNA are tran- 
scription, which creates mRNA molecules at a certain rate, 
and mRNA degradation, which destroys molecules at a 
rate proportional to their abundance. This dynamics gives 
rise to the well-known linear equation of motion [13,14] for 
the concentration x of mRNA of a given gene 

d t x = -j]x + f , (1) 

where r/ is a mRNA-specific degradation constant and / 
is the transcription rate. 

To model changes in the transcription rate over time, 
this deterministic dynamics is generalized to a stochastic 
process described by [15-17] 

dtx ee x{t ± A ; } - x{t) = -r, x + / + y/Dj^m ■ (2) 

At this point, the term £(f) describes lack of knowledge 
regarding regulation by gene-specific transcription factors 
whose concentrations vary over time, and possible exper- 
imental noise. We model this term by a random vari- 
able, whose mean and variance are zero and one, respec- 
tively, without loss of generality. Regulation will be in- 
corporated below. The amplitude y/D/A t of this noise 
term is characterized by a 'diffusion constant' D and the 
interval A t between two successive measurements. Equa- 
tion (2) is well-known as the Langevin equation describing 
an Ornstein-Uhlenbeck process [18]. 

The soundness of (2) as a model of mRNA dynamics 
can be probed by separating deterministic and stochastic 
contributions to the dynamics. The drift term ((xt+i — 
Xt)/A t ) is computed by considering all instances of con- 
centration xt = x in a time course. We use the time 
course data by Spcllman ct al. [19], where expression lev- 
els of all yeast genes were measured at intervals ranging 
from 7 to 20 minutes, see Appendix. The expression ra- 
tios (fluorescence measured relative to a reference sample, 
also termed expression levels) were used as measures of 
mRNA concentration. This is valid provided the fluores- 
cence/concentration response is in the linear regime [20]. 
The result for a specific mRNA shown in Fig. 2(a) is com- 
patible with the linear relationship (dtx) — —rjx + /. 
Genome-wide results will be discussed below. Similarly, 
the distribution of the noise £(i) can be deduced from 
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the data, This distribution is shown for a single gene in 
Fig. 2(b) and is found to be compatible with a Gaussian 
distribution. A Kolmogorov-Smirnov test shows that this 
holds also for all other genes. 
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Fig. 2: Statistics of mRNA concentration dynamics: 
Drift and diffusion, a) The average rate of expression level 
change, or drift {(xt+i — xt)/A t ) for a single mRNA species 
is plotted against the expression levels Xt, by binning a time 
course of 35 measurements into 5 bins. The mRNA chosen for 
this example is encoded by gene YN40, which probably codes 
for a cell wall protein [21]. b) Solving (2) for £(t), the statis- 
tics of the stochastic term can be estimated, shown here for 
the example gene YN40. c) To allow comparison of genes with 
different decay constants and transcription rates, rescaled ex- 
pression values q were used to compute the drift of q for all 
genes in yeast, see text. The genome- wide result agrees well 
with the linear decay law (dtq/rj) — —q indicated by the solid 
line, d) The genome- wide distribution of the noise term P(£), 
see text, e) The standard deviation (J^(q) of the noise term is 
shown against the rescaled expression level q. f) The distribu- 
tion P(£t,£t+A) of noise terms at subsequent time intervals t 
and t + A, see text. 

For a Gaussian noise statistics, the equation of mo- 
tion (2) gives the conditional probability of going from 
concentration x t to Xt+i in a given short time step as 



% /27rL>A t 



exp 



At 
2D 



(d t x + rjx t - fY 



Assuming uncorrelated noise, the propagator (3) gives the 
probability for a complete time course of mRNA expres- 
sion levels, x = (x±, x%, . . . , xt), starting from some initial 
point xi, as 



T-l 



7\/,d(x) = G v j tD (xt+i\x t ) 



(4) 



In order to determine the model parameters we now ask 
which parameter values 77, /, D give the highest likelihood 
(4) to a given time course x. This leads to the straight- 
forward maximum likelihood estimates [22] 



(rj*, /*,£>*) = argmax^ ^-p^/^x) 



(5) 



G v j,D{xt+i\xt) 



(3) 



The time intervals At between adjacent time points fre- 
quently differ from one another. This makes parameters 
estimates from (3) - (5) distinct from fitting data to a 
deterministic model with a quadratic penalty of residu- 
als [23-25]. 

The experimental data [19] gives expression level time 
courses for nearly all ~ 6000 yeast genes. This allows to 
probe the model (2) on a genome-wide scale with high 
statistical accuracy. In order to compare the dynamics of 
genes with different values of the parameters rj, /, D we in- 
troduce rescaled expression levels q = ^2/ (Drj)(rjx— f) so 
the Langevin equation (2) reads d t q = —rjq + ^/27//A t £(i) 
and distribution of q in equilibrium is independent of the 
parameters, P(q) ~ exp{— q 2 /2}. 

The average rate of change (drift) of rescaled expres- 
sion levels q between two successive time steps is shown in 
Fig. 2(c) for all genes. The linear relationship (dtq/rj) = 
— q is obeyed quite accurately. Similarly, the distribution 
of the noise term £(t) across all genes, shown in Fig. 2(d), 
rather closely follows a Gaussian distribution. 

At the genome-wide level a further assumption of the 
Langevin equation (2) can be examined. There the noise 
amplitude was taken to be independent of the expression 
level x. To examine this assumption, the standard de- 
viation a^(q) of the noise term is computed at different 
rescaled expression levels q. Fig. 2(e) shows only an in- 
crease in noise variance of 15% over the whole range of 
expression level changes. This increase of the noise am- 
plitude may be due to increased experimental noise at 
higher expression levels, possibly arising from an intensity- 
dependent sensitivity of the fluorescence measurement of 
the microarrays. An alternative scenario for a changing 
noise-amplitude are stochastic fluctuations of the degra- 
dation rate rj. However, such fluctuations would lead to 
a linear multiplicative noise term in the Langevin equa- 
tion (2), that is, a noise amplitude which linearly increases 
with q. The saturation of cr^(q) at high and low values of 
q which is seen in Fig. 2(e) does not support such a lin- 
ear multiplicative scenario. This result - a small, nonlin- 
ear change in the noise amplitude - contradicts stochastic 
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models with purely multiplicative noise which have been 
used in the literature [26-28]. 

A last assumption behind the Langevin equation (2) 
and the corresponding propagator (3) is the statistical in- 
dependence of noise terms at successive time steps. The 
distribution of successive noise terms P (£t+A, £t) across all 
genes is shown in Fig. 2(f). For small values of £, successive 
noise terms are independent, leading to circular contours 
in Fig. 2(f). For very short intervals A between measure- 
ments, not realised in currently available data, one should 
expect to find correlations between the noise terms £(t) 
and £(t + A). Hence the assumption of uncorrelated noise 
will not hold for arbitrarily small sampling intervals. For 
this reason, the short term propagator (3) was used here, 
rather than the full propagator of an Ornstein-Uhlenbeck 
process with delta-correlated noise. In practice, this choice 
makes little difference for the results obtained. 

These results show that the equation of motion (2) with 
an independent Gaussian distribution of noise of constant 
amplitude is a good description of the observed expression 
level dynamics. Some deviations occur at high and low 
expression levels. The origin of these deviations might, 
however, lie in the details of how the experiments [19] 
were taken. One particular source may be non-linearities 
in the relationship between expression levels and mRNA 
concentrations [20]. Future experiments, taken with dif- 
ferent experimental technologies and at shorter sampling 
times, will allow further tests and possible modifications 
of the model. 

mRNA degradation rates. — We now discuss the 
biophysical parameters inferred from the expression level 
time courses of different genes. The decay times t = vT 1 
inferred using (5) for different genes range from 350s to 
9900s with mean 1400s and standard deviation 700s. De- 
cay times are found to differ somewhat across genes with 
different function, with mRNA of genes involved in tran- 
scription regulation being degraded faster than average 
(r = 1200±100s), and those of genes involved in metabolic 
processes being degraded more slowly (r = 1500 ± 50s). A 
simple possible explanation for these different decay times 
is the need to quickly change concentrations of regulatory 
proteins in response to changing external conditions. 

These results agree with experiments [29] , where mRNA 
decay times were measured after halting mRNA produc- 
tion through a heat shock. (This unphysiological condition 
may well lead to different decay rates than under stan- 
dard conditions.) The measured decay times range from 
r = 160s to t = 7800s with mean 1700s and standard de- 
viation 560s. The experimental study [29] also notes the 
relatively slow degradation of mRNA of genes involved in 
metabolic processes, and the faster degradation of those 
involved in regulatory systems. 

Regulatory interactions. — The transcription of 
genes is controlled by the presence of transcription factor 
proteins. Accordingly, the transcription rate / of a target 
gene depends systematically on the concentration of its 



transcription factors. The dynamic equations of different 
genes are thus coupled by regulatory interactions. The 
resulting correlation of expression levels of different genes 
serves as the basis for exploring regulatory interactions. 

We start by considering the effects of a single transcrip- 
tion factor on the mRNA concentration a; of a given gene. 
The key questions are (i) how does the transcription rate 
of a target gene depend on the concentration of its tran- 
scription factor and (ii) can one identify the targets of a 
given transcription factor from their effect on the tran- 
scription rate ? We generalize the stochastic dynamics (2) 
to 

dtx = -r)x + f(y) + /dTA^M , (6) 

with the transcription rate f(y) depending on the concen- 
tration y of the transcription factor at time t. In the fol- 
lowing, we will neglect post-transcriptional regulation and 
take the mRNA expression level of a transcription factor 
as a proxy for its protein concentration. The functional 
form of f(y) for a given mRNA and a given transcription 
factor can be inferred by discretizing the continuous values 
of y into a small number of bins containing an equal num- 
ber of instances, and inferring the corresponding values 
of / at the centres of these bins. An example, involving 
the transcription factor Swi4 and one of its targets, gene 
YN40, is shown in Fig. 3(a). Similar results are found for 
other transcription factors and their targets. 




(a) 



Fig. 3: Transcription regulation by transcription fac- 
tors, a) The functional form of the transcription rate f(y) 
determined by discretizing the transcription factor expression 
levels y into four bins with equal number of incidences and in- 
ferring the corresponding four values of the transcription rate 
/ (solid squares) along with the parameters rj and D. The re- 
sult, shown here for transcription factor Swi4 and target YN40, 
agrees well with the Michaelis-Menten law (8) displayed by the 
solid line. Error bars indicate the point where the likelihood (4) 
has decayed to e -1 ' 2 of its maximum, corresponding approxi- 
mately to one standard deviation, b) The distribution of bind- 
ing energies (relative to their free energies in solution) of the 
Swi4-DNA and the Swi4-Swi6 interactions, see text. 

Fig. 3(a) displays a non-linear relationship between 
transcription factor concentration and the target tran- 
scription rate showing saturation at high concentrations. 
This is expected since at high transcription factor concen- 
tration most binding sites are occupied by a transcription 
factor molecule, and further increases of concentration no 
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longer have a strong effect. Assuming that a transcrip- 
tion factor may either be bound at a binding site, or be 
suspended in solution or bound non-specifically elsewhere 
on DNA, the probability of a given site being occupied is 
given by elementary statistical mechanics [30] as 



P(y) = 



ye 



-0E 



ye 



-f3E 



-8F 



(7) 



and depends on the transcription factor concentration y, 
binding energy E and a free energy F of the transcription 
factor in solution or bound elsewhere. 

Assuming the transcription rate to depend linearly on 
the probability that the binding site is occupied at a given 
time one obtains 



f(y)=fa + 6P(y) , 



(8) 



where /o is a basal transcription rate in the absence of 
transcription factors and d quantifies the change of the 
transcription rate due to transcription factor binding. 
Positive values of S describe enhancers and negative values 
of 5 indicate repressors. The functional form (8) is the cel- 
ebrated Michaelis-Menten kinetics, first studied in the con- 
text of enzymatic reactions nearly a century ago [31] and 
used widely in transcription modelling [14]. Like many 
other targets of transcription factor Swi4, the relationship 
between transcription rate of target YN40 and transcrip- 
tion factor concentration fits the Michaelis-Menten model 
well, see Fig. 3(a). 

One can turn this result on its head and use the 
Michaelis-Menten model (8) and equations (5) - (6) to 
infer the binding parameters (3(E — F),fo,S of a tran- 
scription factor for different genes. We expect that genes 
with a high response Af = /(y max ) - f(y m m) to changes 
in concentration of a transcription factor are targets of 
that transcription factor. We test this relationship for 
the transcription factor Swi4, and consider genes with dif- 
ferent number of Swi4 binding sequences in their regula- 
tory region, see Figure 4. One finds that Af tends to 
increase with the number of binding sequences, indicat- 
ing that Swi4 acts an enhancer of gene expression. The 
10 genes with the highest Af are CDC9, RNR1, YG3N, 
CRH1, YIOl, RAD27, PRY2, CSI2, PMS5, and CDC21. 
With a single exception, all of these genes have at least 
one copy of the Swi4 binding sequence in their regulatory 
region. Furthermore, 8 of the 10 predicted targets have 
been previously found experimentally [32]. 

Many genes are regulated by more than one transcrip- 
tion factor. Higher organisms in particular owe much of 
their complexity to the combinatorial logic implemented 
by the cooperative binding of several transcription fac- 
tors [33] . For instance, the transcription factor Swi4 forms 
a dimer with transcription factor Swi6, with Swi4 binding 
to regulatory DNA and Swi6 attracting the polymerase 
molecule which starts transcription [34], see Fig. 1(a). For 
this situation, the binding probability (7) generalizes to 

V4V6 



p(Af) 




P(Vi, ye) 



gPiEt-Fi+Ee-Fe) 



V4e 



0(E 6 - 



2/4 Ve 



(9) 



Fig. 4: Genes with a higher number of Swi4 bind- 
ing sites respond more strongly to changes in the 
Swi4 expression level. The distribution of the range A f = 
|/(y ma x) — /(ymin)| of the force term is shown for different sets 
of genes: all genes (black), genes with one or more Swi4 binding 
motifs (red), genes with two or more motifs (green), genes with 
three or more Swi4 binding motifs (blue). Genes with a higher 
number of binding sites are seen to have a larger amplitude of 
the driving force. 



where j/4 and ye are the concentrations of Swi4 and Swi6, 
respectively. The constant E4 gives the energy of Swi4- 
DNA binding and F4 the free energy of Swi4 in solution, 
Eq is the energy of the Swi4-Swi6 interaction and Fq the 
free energy of Swi6 in solution. The Swi4-DNA energy 
depends on the binding sequence in the regulatory region 
of a given gene and may thus vary from gene to gene. On 
the other hand, the Swi4-Swi6 interaction does not de- 
pend on the target gene. Fig. 3(b) shows the distribution 
of [3{Ei — F4) and {3(Eq — F$) inferred by maximum like- 
lihood for different targets of Swi4~Swi6. As expected, 
the variance of the Swi4-Swi6 interaction is significantly 
smaller than that of the Swi4-DNA interaction. 

Discussion. — Some of the properties of gene expres- 
sion are likely to be universal and occur in many organisms 
under many circumstances. An example is the biophys- 
ical implementation of transcription regulation through 
the binding of molecules, with the resulting saturation of 
bound molecules at high concentrations. Another is the 
independent degradation of mRNA molecules, leading to 
the linear drift of mRNA concentrations. (This is less uni- 
versal, since at very high concentrations mRNA molecules 
may outnumber RNase molecules responsible for mRNA 
degradation.) 

Other properties of gene expression are specific to a 
given organism. Indeed, higher organisms, which share 
largely the same set of genes, derive most of their diversity 
from specific changes in how those genes are regulated [35] . 
These specific properties include the transcription factor 
binding sites in the regulatory region of a gene, and the 
interactions between transcription factors. In this paper 
we have used simple stochastic models of mRNA dynamics 
based on the universal properties to infer some of the spe- 
cific properties of gene expression in yeast. Current and 
future developments in molecular biology will bring dy- 
namically resolved data on other types of molecules apart 
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from messenger RNA: concentrations of regulatory pro- 
teins themselves [25,36] or of /iRNA, which silences spe- 
cific niRNA molecules, leading to new challenges in the 
non-equilibrium dynamics of genetic regulation. 
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Appendix. — The mRNA expression data of Spell- 
man [19] comprises four time courses corresponding to dif- 
ferent ways of synchronizing yeast cell cultures to the same 
phase in the cell cycle. Three timecourses were used here 
(termed alpha, cdcl5, cdc28 in [19]), each set equally con- 
tributing to the likelihood (4) . These different timecourses 
correspond to different ways of synchronizing a population 
of cells in a given phase of the cell cycle. The fourth set 
(clu) was not used as it systematically lead to decreased 
likelihoods in (4). The likelihood (4) also allows to iden- 
tify probes on the chip which give results compatible with 
uncorrelated random data; only probes where (4) gave log- 
likelihoods larger than 1.1 the log-likelihood of the data 
under an uncorrelated Gaussian model were used. Set al- 
pha contains 18 samples spaced 7 minutes apart, cdcl5 24 
samples 10 — 20 minutes apart, cdc28 contains 17 samples 
spaced 10 minutes apart. 

The functional classification of genes followed the Gene 
Ontology [37] using the terms transcription regulation ac- 
tivity (GO:0030528) and metabolic process (GO:0008152). 

The canonical binding sequence for Swi4 is "CR- 
CGAAA" where R stands for either G or A [38]. 
Regions 500 base pairs before the transcription ini- 
tiation point were searched for copies of this se- 
quence (or its reverse, its complement, or its re- 
verse complement) . Sequences were taken from 
www . ncbi . nlm . nih . gov/ CBBresearch/Landsman/Cell_ 
cycle_data/upstream_seq. html 24% of all yeast genes 
contain at least one copy of the Swi4 binding sequence in 
their regulatory region, 4% contain two copies or more. 
The YEASTRACT-database [32] list 170 Swi4 targets 
verified experimentally, or about 3% of the yeast genes. 

The example gene YN40 in Figs. 2 and 3(a) was chosen 
due to the high number of 3 Swi4 binding sites within 500 
basepairs and 5 sites within 1000 basepairs of the tran- 
scription initiation site. 



REFERENCES 

[1] Lee T. I. et al, Science , 298 (2002) 799. 
[2] Harbison C. T. et al, Nature , 431 (2004) 99. 
[3] Hertz G. and Stormo G., Bioinformatics , 15 (1999) 
563. 

[4] Bussemaker H., Li H. and SlGGlA E. D., Proc. Natl. 
Acad. Set. USA , 97 (2000) 10096. 



[5] Chua G., Robinson M. D., Morris Q. and Hughes 

T. R., Curr Opin Microbiol , 7 (2004) 638. 
[6] Gene expression omnibus (GEO) http://www.ncbi.nlm. 

nih.gov/geo/ (2007). 
[7] Margolin A. et al, BMC Bioinformatics , 7 (2006) S7. 
[8] HERTZ J., citeseer.ist.psu.edu/ 

hertz98statistical.html (1998). 
[9] Friedman N., Science , 303 (2004) 799. 
[10] Bar- Joseph, Z., Bioinformatics , 20 (2004) 2493. 
[11] Lebre S., http://arxiv.org/abs/0704.2551vl (2007). 
[12] Chu T., Glymour C, Scheines R. and Spirtes P., 

Bioinformatics , 19 (2003) 1147. 
[13] Monod J., Pappenheimer, Jr A. and Cohen-Bazire 

G., Biochim. Biophys. Acta , 9 (1952) 648. 
[14] Alon U., An Introduction to System Biology: Design 
Principles of Biological Circuits (Chapman & Hall, Boca 
Raton, FL) 2007. 
[15] PAULSSON J., Nature , 427 (2004) 415. 
[16] Ozbudak E., Thattai M., Kurtser I., Grossman A. 
and VAN Oudenaarden A., Nature Genetics , 31 (2002) 
69. 

[17] Chen W., England J. and Shakhnovich E., http:// 
arxiv . org/abs/q-bio . MN/0402021 (2004) . 

[18] VAN Kampen N., Stochastic Processes in Physics and 
Chemistry (Elsevier Science, Amsterdam) 1992. 

[19] Spellman P. T., et al, Mol. Biol. Cell , 9 (1998) 3273. 

[20] Carlon E. and Heim T., Physica A, 362 (2006) 433. 

[21] Kawahata M., Masaki K., Fujii T. and Iefuji H., 
FEMS Yeast Res , 6 (2006) 924. 

[22] Honerkamp J., Stochastic Dynamical Systems: Con- 
cepts, Numerical Methods, Data Analysis (VCH, New 
York) 1994. 

[23] Wahde M. and Hertz J., Biosystems , 55 (2000) 126. 

[24] Ronen M., Rosenberg R., Shraiman B. I. and Alon 
U., Proc Natl Acad Sci U S A , 99 (2002) 10555. 

[25] Khanin R., Vinciotti V. and Wit E., Proc. Natl. Acad. 
Sci. USA , 103 (2006) 18592. 

[26] Chen K., Wang T., Tseng H., Huang C. and Kao C, 
Bioinformatics , 21 (2005) 2883. 

[27] Climescu-Haulica A. and Quirk M. D., BMC Bioin- 
formatics , 8 (2007) S4. 

[28] Stokic D., Hanel R. and Thurner S., http://arxiv. 
org/ abs/071 1.4775, (2007). 

[29] Wang Y., et al, Proc. Natl. Acad. Sci. USA , 99 (2002) 
5860. 

[30] Gerland U., Moroz D. and Hwa T., Proc. Nat. Acad. 

Sci. USA , 99 (2002) 12015. 
[31] Michaelis L. and Menten M., Biochem. Z. , 49 (1913) 

333. 

[32] YEASTRACT http: //www. yeastract.com/ (2007). 
[33] Buchler N., Gerland U. and Hwa T., Proc. Natl. 

Acad. Sci. USA , 100 (2003) 5136. 
[34] Andrews B. J. and Moore L. A., Proc Natl Acad Sci 

U S A ,89 (1992) 11852. 
[35] King M. and Wilson A., Science , 188 (1975) 107. 
[36] Bar-Even A.,et al, Nature Genetics , 38 (2006) 636. 
[37] The Gene Ontology Consortium, Nature Genet. , 25 

(2000) 25. 

[38] Chen G., Hata N. and Zhang M., Nucleic Acids Res. , 
32 (2004) 2362. 



p-6 



